Single-particle versus pair condensation of hard-core bosons with correlated hopping 
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We investigate the consequences of correlated hopping on the ground state properties of hard- 
core bosons on a square lattice as revealed by extensive exact diagonalizations and quantum Monte 
Carlo simulations. While for non interacting hard-core bosons the effective attraction induced by 
the correlated hopping leads to phase separation at low density, we show that a modest nearest- 
neighbor repulsion suppresses phase separation, leading to a remarkable low-density pairing phase 
with no single particle Bose- Einstein condensation but long-range two-particle correlations, signaling 
a condensation of pairs. We also explain why the unusual properties of the pairing phase are a real 
challenge for standard one-worm quantum Monte Carlo simulations. 

PACS numbers: 05.30.Jp, 03.75.Kk, 03.75.Lm, 03.75.Hh 
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I. INTRODUCTION 

Models of interacting bosons arise naturally in differ- 
ent fields of condensed matter and of atomic physics. 
The interplay of kinetics and interaction leads to in- 
teresting physical phases like superfluids (SF), Mott in- 
sulators (MI) or supersolids. Phase separation (PS) is 
also a common feature in bosonic systems. Usually, the 
nearest-neighbor hopping and local repulsions are the 
dominant processes. In frustrated physical systems how- 
ever, correlated hopping can be the dominant source of 
kinetic energy. In a mean-field analysis, it has been 
shown that such purely kinetic terms favor the existence 
of a new pairing phase which shows a molecular con- 
densate ((&f&j) ^ 0) without single-particle condensa- 
tion ((oj) = 0) 1 -. Similar phases have been discussed in 
atomic physics in the context of ultracold atoms 2 -^^. 
In particular, a continuous transition from an atomic to 
a molecular condensate is of interest since it would be 
a realization of an Ising transition between two different 
supcrfluid phases. This quantum phase transition has 
to be contrasted with the smooth BEC-BCS crossover as 
e.g. in an atomic Fermi gas. Here both phases can be dis- 
tinguished by their symmetry. In a standard supcrfluid 
the U(l) symmetry is completely broken while in a paired 
superfluid a residual discrete Z2 symmetry remains. Evi- 
dence beyond mean-field for the existence of such a paired 
superfluid (PSF) is still lacking however, and its detailed 
physical properties and the quantum phase transitions 
out of such a phase remain to be understood. 

In this work we study numerically by exact diagonal- 
ization (ED) and quantum Monte Carlo (QMC) the phys- 
ical properties of this exotic pairing phase of bosons. We 
first give an introduction to the model. In Scct. lIIll we dis- 
cuss the ED results, including the global phase diagram 
of the model, the low-density limit, and the relevant cor- 
relation functions. Next we focus on the determination of 
the physical properties for larger system sizes obtained 
by QMC (Sect. HV|) . At the end we collect the major 
findings and give some conclusions. 



II. MODEL 

We consider a model of hard-core bosons with corre- 
lated hopping and nearest-neighbor repulsion on the two- 
dimensional square lattice 
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where rii = b\b i is the boson density at site i, [i is the 
chemical potential, t is the nearest-neighbor hopping am- 
plitude, t' denotes the strength of the correlated hopping, 
and V is a nearest-neighbor repulsion. In the following 
we set t + 1' = 1 and plot the results as a function of t' 
only. We treat the bosons as hard-core particles, i.e. the 
possible local states are restricted to rii <E {0, 1}. 

For vanishing t' , the model can be mapped onto a XXZ 
model in a magnetic field. The phase diagram of the 
model is rich and rather well understoodi^^. At half 
filling, one has a SF for V < 2t and a checkerboard solid 
for V > 2t. For V = 2t the model maps to the Hciscn- 
bcrg model. Away from half filling, a large negative or 
positive chemical potential leads either to an empty or 
a full system. The transition between the checkerboard 
solid and the SF is likely to be first order except at the 
Heisenberg point at half filling. A supersolid phase is not 
realized for this model^. 

For finite correlated hopping, the situation is much 
less clear. In this work we are not interested in pos- 
sible charge ordering phases resulting from the nearest- 
neighbor interaction. We focus on the role of the cor- 
related hopping for the realization of unconventional SF 
phases. For vanishing repulsion, mean field theory target- 
ing a homogeneous solution predicts a phase transition 
between a standard SF with single-boson condensation 
and an exotic PSF with boson-pair condensation for large 
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enough correlated hopping^. Indeed, the basic physical 
effect of the correlated hopping is to act as a binding 
mechanism for bosons. This however opens the possi- 
bility that the system undergoes PS, similar to bosonic 
ring-exchange models^. We will show in this work that 
this is indeed the case. The role of the nearest-neighbor 
repulsion is therefore crucial. We propose the following 
scenario. The system displays PS if the nearest-neighbor 
repulsion is small. On intermediate scales of the repulsion 
the PS is significantly reduced and the boson-pair con- 
densation region is stabilized. Note that this is different 
for bosonic ring-exchange models where a finite repulsion 
only gives rise to a standard SF— . For large repulsion 
the pairs themselves get broken up and a standard SF is 
recovered. 



III. EXACT DIAGONALIZATION 

In this section we discuss results obtained by ED on 
finite clusters. The phase diagram has been obtained 
on systems with up to N s = 36 sites. The low-density 
limit (up to 4 particles) was calculated up to a cluster 
with 7V S = 58. In the following we will first discuss the 
full phase diagram, then the low-density limit and finally 
relevant correlation functions. 

The phase diagram of the model is shown in Fig. [1] as 
a function of t' and density n for different values of the 
nearest-neighbor repulsion. White regions represent SF 
regions, grey regions PS, and black regions the PSF. The 
different phases have been determined in the following 
way. Depending on the chemical potential, the system is 
filled by a certain number of particles. In a conventional 
phase the number of particles increases by one when the 
chemical potential is varied. In the PSF the insertion of a 
single boson is suppressed due to the presence of a pairing 
gap, and pairs of bosons enter the system. In the case 
of PS, boson-pairs are unstable towards cluster forma- 
tion and the density jumps by a finite amount when the 
chemical potential is changed. The three different phases 
SF, PSF, and PS can then be characterized by the way 
in which the total particle number iV p = nN s changes 
upon changing the chemical potential with A(iV p ) = 1, 
A(N p ) = 2, or A(7Vp) > 2 respectively. Typical curves 
of the total number of particles N p as a function of /j, are 
shown in the insets of Fig. [T] 

Let us look first at V = 0. The phase diagram is shown 
in Fig. QJi for the N s = 25 cluster. For not too large val- 
ues of t' or for large densities a SF phase is realized. The 
SF becomes unstable to PS at lower densities for large 
correlated hopping. The tiny shaded region at the left 
edge of the PS at low densities (also present in Fig. [TJd) 
represents a system with exactly two bosons which are 
paired by the correlated hopping. It does not present a 
true low density phase. Practically, the PSF is absent for 
this case. One therefore finds that the correlated hopping 
binds the bosons so strongly that the system phase sep- 
arates between either an empty system or a system at a 
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FIG. 1: Phase diagram obtained by ED. (a) V = and iV s = 
25 (b) V = 2.0 and iV s = 25 (c) V = 2.8 and iV s = 32. 
White represents the SF phase, grey denotes PS, and black 
represents the PSF. The tiny shaded region marks a region of 
pairing which is only present for a finite size system. Insets: 
Total particle number N p as a function of /j, for t' = 0.95. 



finite density. 

The effect of the nearest-neighbor repulsion is shown in 
Fig. [Tb for V = 2.0 and N s = 25 and Fig. QJ for V = 2.8 
and N s = 32. It competes with the attractive interaction 
originating from the correlated hopping. The idea is to 
destabilize the " many-boson" bound states while keeping 
stable bound states of two bosons. In the limit of large 
V, both PS and PSF are destroyed and a standard SF 
at low density is realized. Note also the appearance of 
yet another region of phase separation around densities 
of n = 0.5 for V = 2.8. This region of PS is related to 
the first order transition out of the checkerboard solid 
reported in Ref. At intermediate repulsions, the situ- 
ation is more complex. In Fig. [TJd for V = 2.0, the PSF is 
stabilized for correlated hopping dominating the normal 
hopping at small densities. Still a region of PS remains 
between the PSF and the SF, leading to a region of ex- 
cluded densities between the two stable phases. At the 
even larger value of V = 2.8, (Fig. QJ) the PS is com- 
pletely suppressed, while a reduced region is maintained 
where the PSF is stable. 

The major outcome of the phase diagram for the finite 
clusters is that the PSF is indeed realized in a small win- 
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FIG. 2: Extension of the different phase in the low-density 
limit. White represents the SF phase, grey denotes PS, and 
black represents the PSF. The tiny shaded region marks a 
region of pairing which is only present for a finite size system. 



dow of finite nearest-neighbor repulsion and at low densi- 
ties. In the following we will study the dilute low-density 
limit in order to see whether we get a consistent phys- 
ical picture. We therefore calculated the ground state 
energies £7(1), £7(2), and £7(4) for a 7V S = 58-^ cluster 
with one, two and four particles, and form the two and 
four particle binding energies A2 = £7(2) — 2£7(1) and 
A 4 = £7(4) - 2£7(2). We attribute a bound state of four 
particles (A4 < 0) to a precursor of PS, a bound state of 
two particles (A2 < 0, A4 > 0) corresponds to the PSF, 
and the case of no bound states (A2 > 0, A4 > 0) cor- 
responds to the SF phase, in analogy to the binding of 
hole-pairs in the 2D t— J models. In Fig. [2] the results 
of such an analysis are displayed as a function of t' and 
V. The results arc along the lines we have deduced from 
the calculations at finite density on the smaller clusters 
above. At zero V there is only a SF phase and PS. But 
the PSF can be stabilized by finite V at large correlated 
hopping. For t = 0, a pair of nearest-neighbor bosons 
hops on an effective square lattice with amplitude t', and 
the pair can gain at most a kinetic energy of 4t' equal 
to half the band-width. So the boundary at V = 4 for 
t' = 1 is an exact result. 

Up to now we have only discussed energies and densi- 
ties, which provided evidence in favor of the presence of a 
region in the phase diagram where the bosons form pairs. 
Yet we have not shown that these pairs actually condense. 
Conceptually, one now has to study the asymptotic be- 
havior of the relevant Green's functions, i.e. two-point 
correlation functions for the single-boson SF and four- 
point correlation functions for the PSF. An exponential 
decay of the single-particle Green's function at zero tem- 
perature signals short-range two-boson correlations and 
therefore the absence of a single-boson SF condensate. 
But how does one deal with the PSF? 

It was realized by Yangii that the reduced density ma- 
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TABLE I: General behavior of the maximum eigenvalues A 
for the different superfluid phases SF and PSF 



(i) 

max 



and A m i x 



trices 



P {1) 

(2) 



(2) 
(3) 



are the relevant quantities to investigate in order to 
answer this question. Here pf) is a N s x A s -matrix 
with i, j G {1, .. ., N s } and p?- a A s 2 x A s 2 -matrix with 
i= (*i,*2) 6 {1,...,AT S 2 } and j - (j u j 2 ) e {1,...,A 2 }. 
The behavior of the largest eigenvalue Amax of p^ and 
Amax of p^ as function of the system size N s for con- 
stant density is decisive. For the SF phase one finds 
Amax oc A s and Amax oc A 2 . In contrast for the PSF one 

oc N s . The general behavior is 



has A m ax OC 1 and Amax 



(1) 



summarized in Tab. [T] The different behaviors of A 

(2) 

and Amix can be understood as follows. Let us consider 
first Pij- For any phase without a usual single-boson 
condensate all large distance two-point correlations are 
generically exponentially suppressed, i.e. all entries in 
p\j with \i — j\ ^> 1 arc very small and the matrix has 
only a finite number of sizable entries for fixed j. In the 
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for V = 2 



FIG. 3: The maximum eigenvalues A„ 
and n = 0.1 as a function of N s for t' = (upper left panel) 
and t' — 1 (upper right panel). 

Lower panel: Exponents ot\ and ai for V — 2 as a function of 
t' for constant density n = 0.1. The exponent ai is shown as 
solid curves and the exponent 02 is shown as dashed curves. 
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limit of no correlations when i ^ j, one has a diagonal 
matrix with maximum eigenvalue of order 1. In contrast, 
for a standard SF one has long-range two-point correla- 
tion functions and all entries of p,-7 are of the same order. 
For the case where all correlations are exactly the same, 
one can easily convince oneself that the maximum eigen- 
value is of the order N s . The behavior of p^ is more 
complex. For the case of no single-boson and no boson- 
pair condensate all four-point correlations are suppressed 
exponentially and one again obtains in complete anal- 
ogy a maximum eigenvalue of order 1. For the standard 
SF all four-point correlation functions are long-ranged 
since the four-point correlators contain always two-point 
correlation functions, i.e. all entries of p^ are of the 

(2) 

same order and the maximum eigenvalue A m ax is of or- 
der iV s 2 . The PSF has suppressed two-point correlation 
functions and simultaneously long-range correlations be- 
tween two boson pairs. Therefore all entries involving 
large distances are basically zero except the ones where 
two distances \i\ — jx\ < £ and \i 2 — J2I < £ (or the other 
combination of indices) are within a pairing length. This 

(2) 

gives a maximum eigenvalue A m ax of order N s . In this 
way one can either have long-range order of single bosons 
or long-range order of boson pairs. It docs not make sense 
to speak about ordering in both channels simultaneously 
since a condensate of single bosons is automatically ac- 
companied by a pair condensate. 

In order to compare the behavior of Amax and Amax 
for the different phases, we concentrate on systems with 
an even number of bosons at constant and low densi- 
ties. Since we are limited by the size of the considered 
systems we use the following procedure to characterize 

Amax 

and Amax- We consider the maximum eigenval- 
ues Amax and Amax as a function of N s and N p where 
we restrict N p to be even. Using a linear interpolation 
between the calculated points we obtain the maximum 
eigenvalues Ami x (^Vs, n = const) and Xmix(N s , n = const) 
as a function of N s for a constant density. The result- 
ing data for V = 2,n = 0.1, and the extreme cases 
t' G {0; 1} are shown in the upper panels of Fig. [31 Then 

we fit the functions Amax (N s ,n = const) with aifiN" 1 ' 3 . 
The exponents ol\ for Xmlx(N s ,n — const) and a 2 for 
Xmlx(N s ,n = const) can then be easily used to charac- 
terize the asymptotic behavior of Amix and Amax- The 
exponents a\ and a 2 are shown in the lower panel of 
Fig. El for V = 2.0 and as a function of t' for a density 
n = 0.1. 

The system realizes a SF phase for small values of t'. 
This is nicely reflected in the displayed exponents. ol\ is 
close to the expected value of 1 and a 2 is close to the value 
of 2. Deviations are presumably due to finite size effects. 
It can be clearly seen that these ratios change drastically 
in the regime of large correlated hopping. The devia- 
tions are strongest for t' w 1. Here one obtains ot\ ~ 
and a.2 ~ 1, i.e. values which are close to the expected 
values for the PSF. In the intermediate /j'-regime where 



one has PS no particular exponents are expected or ob- 
served. The investigation of the maximum eigenvalues of 
and p( 2 ) yields therefore additional evidence for the 
presence of the PSF for large t' at low densities, i.e. the 
system realizes a stable condensate of boson pairs. 



IV. QUANTUM MONTE CARLO 

In this section we exploit the results obtained by QMC 
and we complement the physical picture deduced from 
the ED results, i.e. we will establish the existence of the 
PSF in a finite window for finite nearest-neighbor repul- 
sion. The second part of this section will address the 
question of possible phase transitions out of the PSF. 

The model under study does not suffer from a sign 
problem and is therefore accessible to quantum Monte 
Carlo (QMC) techniques. However, the treatment of 
three-site terms is not standard. We will use a SSE 
algorith m 15 i 16 i 17 to tackle the problem as it has been 
shown to be able to treat multi-site termsi£. Some de- 
tails of the algorithm and its implementation are given 
in the Appendix. 

We start the discussion about the SSE results by im- 
portant comments on the numerical efficiency of the one- 
worm SSE algorithm for the case of the pairing phase. 
A basic property of the PSF are exponentially decaying 
two-point correlation functions. Since typical worm algo- 
rithms are based on the insertion of two discontinuities 
into the system which propagate in an extended configu- 
ration spaced (relevant for measuring two-point correla- 
tion functions), the worm updates in space and time are 
short-ranged. Normally, this does not cause a problem 
because the single-boson SF becomes unstable to insulat- 
ing phases which have short ranged Green's functions like 
charge-ordered states. For the PSF this is different. It 
is a long-range ordered phase with short range two-point 
correlation functions. The very nature of the PSF itself 
therefore leads to intrinsic algorithmic problems for sin- 
gle worm updates. Nevertheless, we believe that a certain 
number of important observations can be made. 

The stiffness of a SF phase is determined in SSE by 
measuring the winding numbe r 16 i 20 . If a macroscopic 
fraction of the particles winds around the system, the 
system acquires a finite stiffness. The SSE algorithm uses 
a starting configuration with fixed zero winding number. 
In the single-boson SF the worm algorithm is efficient 
and the algorithm is able to let a large number of bosons 
to wind around the system. In the PSF however, we 
generically expect pairs of bosons to wind around the 
system. This expectation is indeed confirmed by SSE 
simulations on small systems, where winding number his- 
tograms show a very clear even-odd effect, i.e. only even 
winding numbers occur. But the accumulation of a large 
number of winding bosons pairs is exponentially sup- 
pressed by the short range nature of the single worm up- 
dates in the PSF. The higher winding numbers and there- 
fore the stiffness are not well explored in the PSF, based 
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FIG. 4: (a) Density as a function of chemical potential for 
T = 0.05, t = 0.1, f = 0.9, V = 2.0 obtained by QMC for 
7V S = 100. Error bars are one standard deviation. Number 
of seeds is 16. The dotted black line denote ED results with 
iV s = 36. (b) The ratio G(r = max)/n with r = t/50 as a 
function of the chemical potential for the same parameter set 
giving a measure for the one-particle condensate. The region 
between the two grey lines (guide to the eyes) signals the 
finite-temperature coexistence region (COEX). 

on the present worm update scheme. We are currently 
working on two-worm update schemes which should im- 
prove the efficiency of the algorithm in the PSF—. 

In contrast to the winding number, the density of the 
system can be randomly set in the starting configuration. 
It is therefore possible to use for a given parameter set 
a large number of different seeds and average over the 
different runs. In this way one has access to the total 
density and also to density histograms. The latter counts 
how many times a system with a certain total number of 
bosons is realized. One expects a pronounced even/odd 
effect in the PSF which should be absent for a standard 
SF. 

The third quantity we focus on is the equal-time one- 
particle Green's function G(r,r = 0) = (bj.b ). We ex- 
pect that these correlations are exponentially decaying 
in the PSF, due to the pairing gap, while they are long- 
ranged in the conventional SF. 

Based on the ED results we expect that the PSF is real- 
ized at low densities for large correlated hopping and at 
corresponding values of the nearest-neighbor repulsion. 
In the following, we will concentrate on two parameter 
sets: t = 0.1, t' = 0.9, and V = 2.0 on one hand, and 
t = 0.05, t' = 0.95, and V = 2.8 on the other hand. In 
the first case we expect the PSF at low densities to be 
separated from the standard SF by an excluded density 
region of PS. For the second parameter set ED results 
predict the PSF at low densities with possibly a direct 
transition to the standard SF. In the following we always 
investigate the case T = 0.05 

Since the algorithm is grand canonical, we first take 
a look at the total density as a function of the chemical 



FIG. 5: (a) Density as a function of chemical potential for 
T = 0.05, t = 0.05, t' = 0.95, V = 2.8 obtained by QMC 
for iV s = 100. Error bars are one standard deviation. The 
number of seeds is 32. The dotted black line denote ED results 
with N B = 36. (b) Ratio G(r = max)/n with r = \/50 as a 
function of the chemical potential for the same parameter set 
giving a measure for the one-particle condensate. 



potential for both parameter sets in order to determine 
the low-density region of the system. This is shown in 
Fig. Hi and in Fig. [5k, for a N s = 100 system. The dotted 
lines represent the ED results at zero temperature for the 
N s = 36 cluster. We expect from the ED results the PSF 
to be stable for n < 0.11 (fj, < —0.77) for the first set and 
for n < 0.17 (fj, < —0.27) for the second set of parameters. 
Clearly, both parameter sets show no plateaux in n(/x), 
ruling out the realization of an insulating phase. 

For a given system size, a convenient measure of the 
presence of a one-particle condensate is given by the ra- 
tio of the one-particle Green's function at the maximal 
distance allowed by the cluster to the density G(r — 
max)/n. This ratio is plotted in Fig. 2)d and Fig. 03. 
For large chemical potential this quantity is finite and the 
system is in a standard SF phase. On the contrary, if this 
quantity vanishes it signals the breakdown of the SF. The 
finite region between these two limits marks the finite- 
temperature coexistence region (COEX). Clearly, the SF 
is unstable for both parameter sets at low densities. The 
corresponding one-particle Green's functions G(r) show 
a strong decay with distance of the two-particle corre- 
lation functions in this density regime. But the phase 
present at low densities has to be still characterized since 
it could in principle be cither the PSF or PS. To this end, 
the density histogram is displayed in Fig. [6] and Fig. [71 It 
is expected that this histogram is fundamentally differ- 
ent for the two phases. In the PSF one gaussian with a 
large odd/even effect should be present while one expects 
a superposition of two gaussians (one reflecting a subsys- 
tem realizing a SF phase and a second one reflecting a 
PSF) for a PS phase. Focusing only on low densities 
(small values of the chemical potential) the histograms 
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FIG. 6: Density histogram for iV s = 100, T = 0.05, t = 
0.1, t' — 0.9, V — 2.0 and different values for the chemical 
potential. Error bars are one standard deviation. 
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FIG. 7: Density histogram for iV s = 100, T = 0.05, t = 
0.05, t' — 0.95, V — 2.8 and different values for the chemical 
potential. Error bars are one standard deviation. 



for both parameter sets display one gaussian with almost 
only contributions with an even number of particles. This 
establishes that at low densities the PSF is a stable phase. 

Let us briefly summarize our characterization of the 
large t', low density phase. Whenever we have only a fluc- 
tuating even number of bosons in the system and when 
we observe simultaneously exponentially decaying one- 
particle correlations, it is natural to conjecture that we 
are in presence of a paired superfluid phase. Of course, 
on the basis of these Quantum Monte Carlo results alone, 
we cannot exclude the possibility of a paired phase with 
vanishing stiffness since we presently have no direct ac- 
cess to this quantity, but this seems unlikely on the basis 



of the results of the previous section: The size scaling 

(2) 

Amax oc 7V S strongly suggests pair off-diagonal long-range 
order. Besides, an insulating behaviour can be excluded 
since it would be easily detected as plateaux in the den- 
sity as a function of the chemical potential, which are 
clearly not present in the Quantum Monte Carlo data. 
So altogether the existence of a PSF seems to be on solid 
grounds. 

A second important issue is the nature of possible 
quantum phase transitions out of such an unconven- 
tional PSF. We remind the reader that the two superfluid 
phases (SF and PSF) have distinct symmetries enabling 
the possibility of a true quantum phase transition be- 
tween them. So the transition between the PSF and the 
SF can be either a first order transition, or a second or- 
der transition in the Ising universality classi2^£. The 
available data are not sufficient to answer this question 
in a rigorous fashion. This is mainly due to the above- 
mentioned intrinsic algorithmic problems of the worm 
updates. Besides, to decide between a weakly first or- 
der and a second order transition is numerically hard. 
Nevertheless, important trends can be deduced. 

Let us first discuss the case t = 0.1, t' = 0.9, and 
V = 2.0. Results for this parameter set are shown in 
Fig. [Hand Fig. [5] ED predicts a first order transition be- 
tween the two superfluid phases. A first order transition 
at zero temperature corresponds to a jump of the total 
density as a function of the chemical potential. The fi- 
nite temperature nature of the SSE results will smear out 
this jump slightly. Indeed, the SSE data in Fig. [4] display 
a rather sharp dependence of the density as function of 
the chemical potential close to fj, = —0.77 as expected 
from the ED data. Additionally, the density histogram 
in Fig. [6] displays a finite-temperature coexistence region 
(COEX) of the typical behavior of both superfluids near 
fi = —0.77: There is a gaussian with only even particle 
number at lower densities characteristic for the PSF and 
a second gaussian with even and odd particle number at 
higher density representing a standard single-boson SF. 
It is therefore likely that the transition for this parameter 
set is first order. 

ED predicts a continuous transition for the second pa- 
rameter set t = 0.05, t' = 0.95, and V = 2.8. The 
corresponding results are shown in Fig.[?]and Fig. [5] In- 
deed, the density as a function of the chemical potential 
displays a much smoother behavior as shown in Fig. [5] 
A clear jump is not visible but one can identify an inflec- 
tion point close to fi — —0.33 signaling a change of phase. 
The same trend can be seen in the density histogram in 
Fig. [71 There is a rather smooth evolution of the density 
histogram from the PSF to the SF. But the numerical 
data are not sufficient to decide whether a weakly first 
order (existence of a tiny coexistence region) or a second 
order transition between the two superfluids is present. 
It is likely that the second order transition is only real- 
ized (if at all) in a narrow window for large enough values 
of the nearest-neighbor repulsion. 
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CONCLUSION 



APPENDIX: SSE IMPLEMENTATION 



In this work we studied a model of hard-core bosons 
on the square lattice with correlated hopping. We used 
ED and SSE simulations to establish the existence of 
an unconventional phase of condensed boson pairs at 
zero temperature. In the PSF two-point correlations are 
short ranged while there is long-range order in the four- 
point correlations. In other words the bosons form stable 
molecular pairs which condense. We found that a finite 
nearest-neighbor repulsion is crucial for such a phase. 
The correlated hopping acts as an attractive force be- 
tween the bosons. As a consequence, the system is un- 
stable against PS for small nearest-neighbor repulsion. 
A macroscopic part of the bosons glue together. In the 
opposite limit of large repulsion all bound states between 
the bosons are destroyed and a standard SF is realized. 
But the major outcome of this work is that there is a fi- 
nite window of nearest-neigbour repulsion where the for- 
mation of a condensate of paired bosons is stabilized. The 
nearest-neighbor repulsion is large enough to destroy the 
PS but it is small enough not to destroy the boson pairs. 

A second important aspect is the nature of the phase 
transition out of the PSF. This can be either of first or 
second order. Clearly, the realization of a continuous 
transition would be of great interest since it is expected 
to be in the Ising universality class. The detailed nature 
of the transition at zero but also at finite temperature is 
an open issue and can be studied microscopically in our 
model with correlated hopping. The numerical efficiency 
of the algorithm we have used is currently not sufficient 
to answer the question of the nature of the transition 
rigorously. But it is likely that the continuous transition 
is only realized at large values of the nearest-neighbor 
repulsion. When the repulsion is too small, the transition 
is probably first order. 

In conclusion, we have established in this work the 
existence of an unconventional pairing phase of hard- 
core bosons in a microscopic model with correlated hop- 
ping. We have seen that current SSE simulations based 
on worm updates are not efficient in a phase of boson- 
pair condensation. The development of updates schemes 
bases on two worms are currently under study2i and 
should enable the clarification of the nature of the phase 
transitions out of this exotic phase at zero and at finite 
temperature. 



In this appendix we mention some technical points al- 
lowing the reader to follow the numerical implementation 
of the problem. For some introductory literature on SSE 
we refer to Re h 24 ' 25 . The SSE simulations were done 
using a modified version of the ALPS SSE code^ 2 - based 
on the ALPS libraries 2 ^. Conventionally, SSE algorithms 
use bond terms as basic building blocks. The presence 
of three-site terms demands a larger building block. We 
have chosen as basic objects of our problem operators 
acting on a plaquette of the square lattice^. In the fol- 
lowing we denote the four sites of one plaquette i, j, k, 
and / in anti-clockwise direction putting the first site in 
the left lower corner. We have three classes of vertices. 
First, there are diagonal vertices 

#l,a = Clijkl (A.l) 

A* 

+ 4 [niljki + rijliki + rikhji + nihjki 
V 

-— [niUjlki + ii, in I,, + n k nil vi + ninJkj] . 

Here C is a constant which has to be adjusted such that 
all diagonal contributions are not negative, i.e. that there 
is no sign problem. Second, there are bond vertices 



(A.2) 
(A.3) 
(A.4) 
(A.5) 



where By = b\b l + b\b r Third, we have vertices origi- 
nating from the correlated hopping 



#2, a = 


2 BijIlk 


H3, a = 


2 B jkhi 


#4, a = 




#5, a = 





^6,a 
^7,a 

H9, a 



= t'Tujlk 

= t'lkUiBij 

= t'Tjikli 

= I'h 11, IS, 1. 

= t'Tkijh 

= t'liUkBij 

= tTuklj 

= t'ljmBik 



(A.6) 
(A.7) 
(A.8) 
(A.9) 
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A second important issue is the choice of the updates. 
Beside standard diagonal updates we use directed loop 
updates. The determination of the travel probabilities of 
the two discontinuities in order to fulfill detailed balance 
for the created directed loops is a non-trivial problem. 
A successfull strategy is to minimize the bounce prob- 
abilities having in mind that one wants to create large 
directed loops. Following KciM-, one can write down a 
solution which minimizes the bounce amplitudes. The 
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bounce processes are only non-zero in a parameter regime dominates the other weights, and bounce free solutions 
where one term in the Hamiltonian dominates the other can be obtained^, 
terms. In all other cases, the maximum weight does not 
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